Assessment: Exploratory Analysis of the Factors Influencing Emergency Contraceptive Prescribing in Scotland

Introduction

Preventing unintended pregnancy is critical for individual well-being and public health. Unprotected sex occurs for multiple, complex reasons including sexual assault, reproductive coercion, contraceptive mishaps, and lapses in contraceptive use Cleland et al. (2014) . Emergency contraception prevents pregnancy following unprotected sex or contraceptive failure. In Scotland emergency contraception is available as either a pill or intrauterine device. Oral contraceptives are the most commonly used, making up 90% of emergency contraceptive prescriptions@Emergency23.

Pharmacies in Scotland dispense two types of emergency contraceptive pills (ECP):

• Levonorgestrel (brand name: Levonelle) - Effective within 72 hours (three days) of unprotected sex
• Ulipristal Acetate (brand name: ellaOne) - Effective within 120 hours (five days) of unprotected sex

This analysis explores the factors influencing prescribing trends for ECP in Scotland from 2019 to 2022, using data from the Public Health Scotland Community Prescriptions Database. By investigating ECP patterns over time, geography and socioeconomic status, this report aims to identify disparities in access and inform strategies to improve public health outcomes.

The analysis addresses the following questions:

1. What are the trends in prescribing of ECP between 2019 and 2022?
2. How does the ECP prescription rate vary by Health Board regions in Scotland?
3. Are pharmacies in more deprived areas prescribing more ECP than those in affluent areas?
4. Is there geographical variation in prescribing practices of ECP by Health Board region?

Load necessary libraries:

library(tidyverse)
library(janitor) 
library(gt) 
library(here)
library(lubridate)
library(sf)
library(ggspatial)
library(plotly)

Read in datasets:

# Read in the Health Board names (HB_names). See code for link.
HB_names <- read_csv("https://www.opendata.nhs.scot/dataset/9f942fdb-e59e-44f5-b534-d6e17229cc7b/resource/652ff726-e676-4a20-abda-435b98dd7bdc/download/hb14_hb19.csv") %>% 
  clean_names() # ensures column names are unique, lower case and spaces and special characters are replaced with underscores. I will use this function for all read ins to ensure my data is consistent and easy to manipulate.

# Read in population data per Health Board from the 2022 census data. Available from: https://statistics.ukdataservice.ac.uk/dataset/scotland-s-census-2022-uv102a-age-by-sex/resource/b2d295c2-af53-4b3d-a075-7815cadd9060 
all_population_data <- read_csv(here("data", "UV103_age_health_board_census.csv"), skip = 10) %>% # locates the csv and excludes the first 10 lines in the csv as they are redundant
  rename(Spare = "...6", # remove unused columns
         hb_name = "Health Board Area 2019",
         hb_population = Count) %>% # rename() formats the data to match the prescriptions dataframe
  filter(Sex == "Female") %>% # filter female population, as men do not take ECP
  select(hb_name, Age, hb_population) %>% # select columns of interest
  mutate(hb_name = paste("NHS", hb_name)) %>%  # change hb_name column format to match the Health Board dataframe format
  clean_names()

# filter dataset for the total female population for each health board:
population_data <- all_population_data %>% 
  filter(age == "All people") %>% 
  select(hb_name, hb_population)

# select the population aged 16 to 34 years as this is the population the population most likely engaging with risk taking behaviors such as unprotected sex. I chose 16 as this is the age of consent in Scotland, and 34 as my cut off for young adults. 
young_population_data <- all_population_data %>% 
  filter(age%in% c(16:34)) %>% # filter ages of interest
  mutate(age = as.numeric(age)) %>% # make numerical so they can be placed in to buckets
  mutate(age_group = case_when(between(age, 16,34)~'16-34')) %>% # make bucket 16-34 years
  group_by(hb_name, age_group) %>% 
  summarise(pop_hb_16_34 = sum(hb_population)) # sum total population aged 16 to 34 per health board region. 

Define a function to read in the prescription data for a defined year:

#For efficiency I created a function to read in the prescriptions datasets. I downloaded 12 months of prescription data for each year from 2019 to 2022. I placed the 12 months of data into their relevant folder named all_months_year, where year was specific to the data it contained.
read_all_prescriptions <- function(year){
  all_files <- list.files(here("data", paste0("all_months_", year)), pattern = "csv") #list.files() retrieves files from the relevant directory, and the paste0() dynamically constructs the folder name based on the year variable placed into the function.
  all_prescriptions <- all_files %>%
    map_dfr(~read_csv(here("data", paste0("all_months_", year),.))) %>% #map_dfr() row-binds the datasets being read in
    clean_names() %>%
    drop_na(bnf_item_description) # drop the rows with missing bnf_item values
  return(all_prescriptions)}

Read in datasets using read_all_prescriptions() function and start data wrangling:

# all_prescriptions_2019 <- read_all_prescriptions(2019)
# all_prescriptions_2020 <- read_all_prescriptions(2020)
# all_prescriptions_2021 <- read_all_prescriptions(2021)
# all_prescriptions_2022 <- read_all_prescriptions(2022)
# 
# # 2019 dataset has hbt2014 as a column name instead of hbt. Rename to make column name consistent
# all_prescriptions_2019 <- all_prescriptions_2019 %>%
#   rename(hbt = "hbt2014")
# 
# # combine 4 years of data into one dataset to make it easier to wrangle
# combined_prescriptions <- bind_rows(
#   all_prescriptions_2019,
#   all_prescriptions_2020,
#   all_prescriptions_2021,
#   all_prescriptions_2022 )%>% # I mutate and select the prescriptions of interest to prevent a very large dataset from being stored in my environment (prevents R slowing down and crashing)
#   mutate(
#     date = parse_date_time(paid_date_month, "ym"), # use lubridate to format date
#     drug_simple = case_when(
#       str_detect(bnf_item_description, "LEVONO") ~ "Levonorgestrel",
#       str_detect(bnf_item_description, "ULIPR") ~ "Ulipristal Acetate",
#       TRUE ~ "Other")) %>% # used case_when() to select ECP even when they have different dosages
#   filter(drug_simple != "Other",!is.na(date)) %>% #remove prescriptions of no interest and any missing date values
#   filter(hbt != "SB0806") %>% # filter out SB0806 as it is not a health board (Scottish Ambulance Service)
#   filter(!is.na(hbt))

# save combined_prescriptions dataset to a csv in my directory, to prevent many large datasets being stored in my environment and slowing down R. 
# write_csv(combined_prescriptions,"data/combined_prescriptions.csv")
combined_prescriptions <- read_csv(here("data","combined_prescriptions.csv"))

Join and wrangle data naming the output variable ECP_scripts:

# join the prescriptions dataset to the Health Board names and health board population datasets
ECP_scripts <- combined_prescriptions %>% 
  full_join(HB_names, by = c("hbt" = "hb")) %>% # join with Health Board names
  full_join(population_data, by = "hb_name") %>% # join with population data
  select(gp_practice, date, drug = drug_simple, hb_name, paid_quantity,hb_population) %>%  # select columns of interest
  mutate(month = factor(month(date), levels = 1:12, labels = month.abb), .after = date) %>%   
  mutate(year = factor(year(date)), .after = month) # extract month and year as a factor with labels to help when plotting

Key Results

What are the trends in prescribing of ECP between 2019 and 2022?

# calculate total population of Scotland
total_female_population <- sum(population_data$hb_population, na.rm=TRUE) 

# calculate the number of prescriptions of each drug per month
ECP_by_month <- ECP_scripts %>%
  group_by(date, drug, month, year) %>%
  summarise(total_quantity_month = sum(paid_quantity, na.rm = TRUE)) %>%
  ungroup() %>% 
  drop_na(drug) %>% 
  clean_names()

# calculate prescription rate per 100,000 people for each drug and month
ECP_trend_plot_data <- ECP_by_month %>%
  group_by(drug, year, month) %>% 
  summarise(prescriptions_per_100000 = (total_quantity_month / total_female_population)*100000)

ECP_trend_plot <- ECP_trend_plot_data %>%
  ggplot(aes(x = month, y = prescriptions_per_100000, group = year, linetype = year, color = year)) + # plot data for each year on same graph
  geom_line(aes(size = ifelse(year == "2020", 0.9, 0.7))) + # make the 2020 line slightly thicker to make it stand out 
  facet_wrap(~drug, scales = "free_y") + # separate the two types of ECP
  scale_linetype_manual(values = c("2019" = "dotdash", "2020" = "solid", "2021" = "dashed", "2022" = "dotted"),name = "Year") + # make each year identifiable by linte type
  scale_color_manual(values = c("2019" = "grey70", "2020" = "#4575b4", "2021" = "grey70", "2022" = "grey70"),name = "Year") + # make Covid-19 year blue (2020 had the longest lockdown, March till May)
  scale_size_identity() + # ensures ggplot uses the size values as outlined in geom_line()
  scale_y_continuous(limits = c(0, NA), breaks = scales::pretty_breaks(n = 5))+ # ensure y-axis starts at 0 and has evenly spaced breaks 
  labs(
    title = "Trends in Emergency Contraception Prescribing 2019–2022 \n in Scotland.",
    subtitle = "Four years of population adjusted prescribing rates\n  highlighting the impact of Covid-19 (blue)",
    x = "Month",
    y = "Prescriptions per 100,000 women") +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 18, hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(size = 14, hjust = 0.5),
    axis.title.x = element_text(size = 12),
    axis.title.y = element_text(size = 12),
    axis.text.x = element_text(angle = 45, hjust = 1, size = 9),
    axis.text.y = element_text(size = 9),
    strip.text = element_text(size = 12),
    legend.title = element_text(size = 12),
    legend.text = element_text(size = 12),
    legend.position = "bottom",
    plot.margin = margin(10, 30, 10, 30)) # increase size of margins

ECP_trend_plot

Figure 1: Prescribing trends of Emergency Contraception from 2019 to 2022

Figure 1 explores the monthly variation in prescribing rate of ECP from 2019 to 2022. 4 consecutive years of data were analysed to investigate whether specific events, such as the Covid-19 pandemic or the Common Wealth Games influenced ECP prescribing rates.

Interesting trends to note include:

• Levonorgestrel is more commonly prescribed than Ulipristal Acetate. Each month an average of 150 prescriptions of Levonorgestrel and 3 prescriptions of Ulipristal Acetate are dispensed per 100,000 women across Scotland.

• The prescribing rate of Levonorgestrel has decreased annually between 2019 and 2022. This trend may reflect public health efforts to promote long-acting reversible contraception (LARC), thereby reducing the need for ECP. (Long?) Alternatively it could result from challenges in accessing ECP during the Covid-19 pandemic.

• A significnat reduction in Ulipristal Acetate prescriptions is observed from 2020 onwards. Due to low prescribing rates of Ulipristal Acetate there is some stochastic variation however the downward trend in 2020 is clear. This likely correspounds to the European Medicines Agency’s 2020 guidance recommending teh withdrawal of Ulipristal Acetate for the treatment for uterine fibroids due to safety concerns. From 2020, prescriptions of Ulipristal Acetate appear specific to ECP use only. (Ulipristal20?)

• Seasonal peaks are obserevd in the summer months of all years. This suggests an increased demand for ECP during what could be described as a ‘Summer of Love’.

• A marked decline in Levonorgestrel prescriptions ocurred between February and May of 2020. This coincides with the Phase 1 Covid-19 lockdown, suggestring that national lockdown may have reduced opportunities for unprotected sexual intercourse.

How does the ECP prescription rate vary by Health Board regions in Scotland?

Wrangle data to get the columns and values of interest:

ECP_anual_rate_data <- ECP_scripts %>% 
  group_by(hb_name, drug) %>% # aggregate data by Health Board and drug
  summarise(total_quantity_4_years = sum(paid_quantity, na.rm = TRUE), # total prescriptions over 4 years
    avg_annual_total_quantity = total_quantity_4_years / 4, # average annual prescriptions
    hb_population = first(hb_population), # hb_population is consistent within each hb_name
    .groups = "drop") %>% # ungroup data
  drop_na(drug) %>% # remove rows with no drug values
  mutate(avg_annual_presc_100000 = avg_annual_total_quantity * 100000 / hb_population) %>% # average annual rate of prescriptions per 100,000
  select(hb_name, drug, avg_annual_presc_100000, hb_population) %>% # select columns of interest before pivot
  pivot_wider(names_from = drug, values_from = avg_annual_presc_100000, names_glue = "{drug}_rate") %>% # Rename columns for clarity
  clean_names() # clean column names following pivot

# calculate total ECP prescription rate per health board
ECP_anual_rate_data <- ECP_anual_rate_data %>% 
  mutate(total_ECP_rate = rowSums(select(., levonorgestrel_rate, ulipristal_acetate_rate), na.rm = TRUE)) %>%  #sum rates 
  arrange(desc(total_ECP_rate)) # arrange by total rate in descending order

#calculate percentage young people per health board to explore if there is an association here 
ECP_anual_rate_data_young <- ECP_anual_rate_data %>% 
  full_join(young_population_data) %>% # to get the population of those aged 16-34 years 
  group_by(hb_name) %>% 
  mutate(prop_young_ppl_hb = pop_hb_16_34/hb_population) %>% # calculate proportion of young people in each health board 
  ungroup()

Plot table 1:

annual_avg_ECP_table <- ECP_anual_rate_data_young %>% 
  select(hb_name, levonorgestrel_rate, ulipristal_acetate_rate, total_ECP_rate,prop_young_ppl_hb) %>% # Select relevant columns
  gt() %>% 
  cols_label(hb_name = "Health Board",
             total_ECP_rate= "Total",
             levonorgestrel_rate= " Levonorgestrel",
             ulipristal_acetate_rate=" Ulipristal Acetate",
             prop_young_ppl_hb = "% Aged 16-34 Years") %>% # rename columns to make reader-friendly
  fmt_number(columns = c(levonorgestrel_rate, ulipristal_acetate_rate, total_ECP_rate, prop_young_ppl_hb), decimals = 0) %>% # no decimal points as false accuracy detracts from the message in the data
  cols_align(align = "center", columns = c(levonorgestrel_rate,ulipristal_acetate_rate,total_ECP_rate, prop_young_ppl_hb)) %>%  # centre column names
  grand_summary_rows(columns = c(levonorgestrel_rate,ulipristal_acetate_rate, total_ECP_rate), fns = list("Overall Average" = ~mean(., na.rm = TRUE)), fmt = list(~ fmt_number(., decimals = 0))) %>% # add an overall average for the rate columns
  fmt_percent(columns = prop_young_ppl_hb, decimals = 0) %>%  # add percentage sign
  tab_header(title = md("**Average Annual Rate of Emergency Contraception Prescriptions by Health Board in Scotland.**"),
             subtitle = md("Rate per 100,000 women, derived from the mean prescription rates across the years 2019 to 2022. Health Boards are ranked in descending order.")) %>% # add a title and subtitle; md() allows text formatting from mark down 
  tab_spanner(label = md("*Prescription rate per 100,000 women*"), columns = c(levonorgestrel_rate,ulipristal_acetate_rate, total_ECP_rate)) %>% # add a title to the rate columns.
  tab_source_note(md("*Data from Public Health Scotland. Available from: (https://www.opendata.nhs.scot/dataset/prescriptions-in-the-community)*")) %>% 
  tab_stubhead(md("**2019-2022**")) %>% 
  tab_footnote(footnote = "includes Capital City, Edinburgh", locations = cells_body(columns = hb_name, rows = 2))%>% 
  opt_stylize(style = 6, color = "blue")

annual_avg_ECP_table

Average Annual Rate of Emergency Contraception Prescriptions by Health Board in Scotland.
Rate per 100,000 women, derived from the mean prescription rates across the years 2019 to 2022. Health Boards are ranked in descending order.
2019-2022	Health Board	Prescription rate per 100,000 women			% Aged 16-34 Years
		Levonorgestrel	Ulipristal Acetate	Total
	NHS Greater Glasgow and Clyde	3,376	29	3,404	26%
	NHS Lothian1	2,495	56	2,551	27%
	NHS Ayrshire and Arran	2,027	15	2,043	19%
	NHS Tayside	1,365	22	1,387	22%
	NHS Forth Valley	1,268	73	1,342	22%
	NHS Grampian	1,210	16	1,226	22%
	NHS Fife	1,049	54	1,103	22%
	NHS Lanarkshire	784	48	832	22%
	NHS Shetland	779	17	796	19%
	NHS Orkney	609	7	616	17%
	NHS Dumfries and Galloway	515	24	539	17%
	NHS Western Isles	192	197	389	16%
	NHS Highland	285	52	337	17%
	NHS Borders	242	76	318	17%
Overall Average	—	1,157	49	1,206	—
Data from Public Health Scotland. Available from: (https://www.opendata.nhs.scot/dataset/prescriptions-in-the-community)
1 includes Capital City, Edinburgh

Table 1 presents the annual prescription rate of ECP per Health Board, alongside the proportion of young people living in each Health Board.

Key points to note:

• NHS Greater Glasgow and Clyde and NHS Lothian have the highest annual ECP prescription rate. These regions encompass Scotland’s largest cities, Glasgow and Edinburgh, which may contribute to higher demand for ECP services.

• Health boards with higher ECP prescription rates generally have a higher proportion of young people. This suggests that young people aged 16 to 34 years are more likely to require emergency contraception, in keeping with current literature.@mehta_etal22

• NHS Ayrshire and Arran is an outlier. Despite having the third highest ECP prescription rate, it has a relatively smaller population of young people (19%). This warrants further exploration, as it may indicate greater use of ECP services by under 16s or over 35s in this region helping to guide targeted interventions.

• Remote and rural Health Boards such as NHS Highland and NHS Borders had the lowest rate of ECP prescribing. They also had a smaller proportion of young people living in their health board, which may explain this observation.

• Levonorgestrel is more commonly prescribed than Ulipristal Acetate across most Health Boards. However, in the NHS Western Isles Health Board the prescribing rates for Levonorgestrel and Ulipristal Acetate are nearly equal This variation may reflect disparities in prescribing practices between Health Boards. This is explored further in figure 3.

Are pharmacies in more deprived areas prescribing more ECP than those in affluent areas?

Read in and wrangle data:

# To measure deprivation I have used the Scottish Index of Multiple Deprivation. 2020v2 dataset from Public Health Scotland. (link available in code)
SIMD <- read_csv("https://www.opendata.nhs.scot/gl/dataset/78d41fa9-1a62-4f7b-9edb-3e8522a93378/resource/acade396-8430-4b34-895a-b3e757fa346e/download/simd2020v2_22062020.csv") %>%
  clean_names() %>%
  select(data_zone, simd2020v2hb_decile)

# I chose to use GP Practices and List Sizes from October 2022, as this was the closest dataset I could find which correlated with the final year of my prescriptions dataset:
gp_addresses <- read_csv("https://www.opendata.nhs.scot/dataset/f23655c3-6e23-4103-a511-a80d998adb90/resource/1a15cb34-fcf9-4d3f-ad63-1ba3e675fbe2/download/practice_contactdetails_oct2022-open-data.csv") %>%
  clean_names() %>%
  select(practice_code, gp_practice_name, data_zone)

# Create ECP_GP by using the GP_addresses dataset to map the GP practice code to a datazone. I then used the column datazone to full_join() the SIMD dataset to the prescriptions dataset.
ECP_GP <- ECP_scripts %>%
  filter(!gp_practice %in% c(99996, 99997, 99998)) %>% # remove dummy GP practice codes as they do not have a known gp practice code so cannot be mapped to a SIMD.
  left_join(gp_addresses, by = c("gp_practice" = "practice_code")) %>%
  #left_join(data_zones, by = "data_zone") %>%
  left_join(SIMD, by = "data_zone") %>%
  drop_na(simd2020v2hb_decile) %>% # need SIMD value to for plot
  group_by(gp_practice) %>% 
  mutate(total_quantity_per_gp = sum(paid_quantity)) %>% 
  clean_names()

# calculate the number of GPs per SIMD (account for each SIMD having a different number of GPs) 
toal_no_GP_per_SIMD <- ECP_GP %>% 
  group_by(simd2020v2hb_decile) %>%
  mutate(unique_gp_count_per_SIMD = n_distinct(gp_practice))

Plot figure 2:

ECP_SIMD_barchart <- toal_no_GP_per_SIMD %>%
  group_by(simd2020v2hb_decile, drug) %>% 
  summarise(prescriptions_gp = (total_quantity_per_gp / unique_gp_count_per_SIMD), # calculate the prescriptions per GP for each SIMD decile and drug
    .groups = "drop") %>%  # Ungroup data after summarisation
  ggplot(aes(x = factor(simd2020v2hb_decile, levels = 1:10),y = prescriptions_gp,  fill = drug, text = paste("Decile:", simd2020v2hb_decile, "<br>Drug:", drug))) + # customise hover box for interactive chart
  geom_col() +
  scale_fill_manual(values = c("#4575b4", "#91bfdb"), name = "Drug Type") +  # add colour palette
  labs(title = "Emergency Contraception Prescriptions by \n Scottish Index of Multiple Deprivation (SIMD) 2019 to 2022, \n normalised by number of GPs per SIMD decile.",
    x = "SIMD Decile \n (1 = Most Deprived, 10 = Least Deprived)",
    y = "Total Prescriptions",
    fill = "Drug") +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(size = 14, face = "bold", hjust = 0.5),
    plot.subtitle = element_text(size = 12, hjust = 0.5),
    axis.title.x = element_text(size = 12),
    axis.title.y = element_text(size = 12),
    axis.text.x = element_text(size = 10, hjust = 1),
    axis.text.y = element_text(size = 10),
    legend.position = "right",
    legend.title = element_text(size = 12, face = "bold"),
    legend.text = element_text(size = 10))

ECP_SIMD_barchart <- ggplotly(ECP_SIMD_barchart, tooltip = "text") # make plot interactive
ECP_SIMD_barchart

Figure 2: Interactive Barchart showing Emergency Contraception Prescriptions by Scottish Index of Multiple Deprivation (SIMD).

Generally as SIMD decile increases the number of ECP prescriptions decreases. GPs in SIMD1, the most deprived areas, prescribe the most ECP, whilst increasing affluence is associated with a decline in ECP prescriptions. This suggests that individuals from more deprived areas utilise ECP services more frequently, highlighting a potential inequality in health needs.

This disparity may be due to factors such as education, awareness of safe sexual practices, environmental influences, or limited access to long-acting reversible contraception (LARC). Further investigation into these factors is essential to address this inequality.

Is there geographical variation in the prescribing practices of ECP by Health Board region?

This analysis explores if there was any correlation between the type of contraception prescribed and the location of the pharmacy. The ratio of total Levonorgestrel prescriptions to total emergency contraception prescriptions was calculated for each Health Board Region using the formula:

\[ \frac{Levonorgestrel}{Levonorgestrel + Ulipristal Acetate} = Ratio \]

A score close to 1 indicates that the majority of prescriptions in this region are for Levonorgestrel, whilst a score closer to 0 suggests that the majority of prescriptions are Ulipristal Acetate. This approach provides a useful way of identifying variation in prescribing practices.

Figure 3:

# load the NHS Health board Shapefile downloaded from learn page
NHS_healthboards <- st_read(here("data", "NHS_healthboards_2019.shp")) %>% 
  mutate(HBName = paste("NHS", HBName)) %>% # format to match ECP_scripts dataset
  clean_names()

## Reading layer `NHS_healthboards_2019' from data source 
##   `C:\Data_science\B273025\data\NHS_healthboards_2019.shp' using driver `ESRI Shapefile'
## Simple feature collection with 14 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 7564.996 ymin: 530635.8 xmax: 468754.8 ymax: 1218625
## Projected CRS: OSGB36 / British National Grid

# calculate the ratio of Lev to Uli prescribed 
variation_ECP_prescribed <- ECP_scripts %>% 
  group_by(hb_name,drug) %>% 
  # calculate the total of Lev and Uli prescribed per health board
  summarise(total_each_drug_type = sum(paid_quantity, na.rm = TRUE)) %>%
  drop_na(drug) %>% 
  #pivot_wider to move drug names to columns 
  pivot_wider(names_from = drug, values_from = total_each_drug_type) %>%
  clean_names() %>% # consistency 
  mutate(levo_to_uli_ratio = levonorgestrel / (levonorgestrel + ulipristal_acetate)) # calculate ratio 

# Join spatial data with variation_ECP_prescribed
variation_ECP_prescribed <- NHS_healthboards %>%
  left_join(variation_ECP_prescribed)

#Create map in ggplot
map_variation_ECP_prescribed <- variation_ECP_prescribed %>%
  ggplot(aes(fill = levo_to_uli_ratio)) +
  geom_sf(size = 0.1, colour = "grey50", alpha =0.9) +
  scale_fill_distiller(palette = "Blues", direction = 1) +
  labs(title = "Geographical variation in Emergency Contraceptive \n Pill prescribing in Scotland.", 
       subtitle = "Heatmap showing Levonorgestrel prescriptions as a proportion of total \n Emergency Contraceptive Pill (ECP) prescriptions by Health Board region", 
       fill = "Levonorgestrel to ECP ratio", caption = "Data Source: Public Health Scotland Prescriptions in the Community, 2019-2022") +
  coord_sf() +
  theme_void() +
  theme(
    plot.title = element_text(face = "bold", size = 18, hjust=0.5), 
    plot.subtitle = element_text(size = 14, hjust=0.5), 
    legend.title = element_text(size = 12),
    legend.text = element_text(size = 9, hjust=0.5),
    legend.direction = "vertical",
    legend.box = "horizontal") +
  annotation_scale(location = "tl") +
  annotation_north_arrow(
    location = "tl",    
    pad_y = unit(0.5, "in"),    
    style = north_arrow_nautical(fill = c("grey40", "white"),line_col = "grey20"))

map_variation_ECP_prescribed

Figure 3: Heatmap exploring geographical variation in Emergency Contraception Pill prescribing

Health Board Regions in central Scotland predominantly prescribe Levonorgestrel, where as more remote and rural Health Board regions, such as the NHS Western Isles, prescribe higher proportions of Ulipristal Acetate. This difference in prescribing practices may reflect disparities in access to emergency contraception in remote and rural regions. Patients in these regions may present later to health services later, prompting prescribers to dispense Ulipristal Acetate due to its longer efficacy window. Alternatively, this variation could be attributed to differences in prescriber preferences or guidelines in more remote areas.

Conclusion

Prescriptions of emergency contraceptive pills (ECP) have decreased over the years 2019 to 2022. In 2020, there was a notable reduction in Levonorgestrel prescriptions, likely influenced by the Covid-19 lockdown. Ulipristal Acetate had a marked decline from 2019 onwards likely reflecting updated prescribing guidance, which discontinued its use in the treatment of uterine fibroids.

Health Board regions with larger populations of young people and major cities had higher rates of ECP prescriptions, whilst more remote regions with a smaller proportion of young people had a lower ECP prescription rate.

ECP prescriptions were more prevalent in areas of greater deprivation.

Geographically, remote pharmacies appeared to prescribe higher proportions of Ulipristal Acetate. This may be due to challenges in accessing pharmacies within 48 hours of unprotected sex, necessitating the use of Ulipristal Acetate, which has a longer efficacy window.

In conclusion, ECP prescribing patterns are influenced by multiple factors including societal events such as the Covid-19 pandemic, age, deprivation, and geography.

Recommendations from analysis

This report highlights significant disparities in the prescribing of emergency contraception across different geographic regions and levels of deprivation in Scotland. Further research is needed to understand the underlying causes of these trends, considering if there are barriers to access or differences in prescriber practices. These insights will be crucial for informing future sexual health policies and interventions aimed at reducing health inequalities and ensuring equitable access to emergency contraception.

Limitations of the dataset

• This report assumes that Levonorgestrel and Ulipristal Acetate are dispensed exclusively for emergency contraception, wihtout considering potential off-licence uses.

• Differences in the volume of ECP prescribed across Health Boards may result from patients living in one Health Board region but obtaining their ECP from a pharmacy in another Health Board region. This may result in inaccuracies in attributing prescription rates to specific Health Boards.

• The prescription dataset included ‘dummy practice codes’ for prescriptions which lacked a specific practice code. This excluded a large number of ECP prescriptions from the ECP-deprivation analysis as it was impossible to link the prescriptions to a datazone and subsequent SIMD decile. This may impact the accuracy of the analysis.

Suggestions for future analysis

• With additional time it would be interesting to design a Directed Acyclic Graphic to visually and mathematically model the probabilistic relationships between variables such as age, deprivation and geography on ECP prescribing patterns in Scotland.

• It would be interesting to conduct statistical analysis to explore the trends in this report.

References

I linked this rmd to Zotero Referencing Software to produce in text citations with an updating bibliography. I chose to format my refernces in BMJ house style as this is commonly used in medical literature.

Cleland, Kelly, Elizabeth G. Raymond, Elizabeth Westley, and James Trussell. 2014. “Emergency Contraception Review: Evidence-Based Recommendations for Clinicians.” Clinical Obstetrics and Gynecology 57 (4): 741–50. https://doi.org/10.1097/GRF.0000000000000056.